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The pseudospectral method is a powerful tool for finding highly precise solutions of Schrodinger's 
equation for few-electron problems. Previously we developed the method to calculate fully correlated 
S-state wave functions for two-electron atoms [J|. Here we extend the method's scope to wave 
functions with non-zero angular momentum and test it on several challenging problems. One group 
of tests involves the determination of the nonrelativistic electric dipole oscillator strength for the 
helium IS — s> 2 l ¥ transition. The result achieved, 0.27616499(27), is comparable to the best in the 
literature. The formally equivalent length, velocity, and acceleration expressions for the oscillator 
strength all yield roughly the same accuracy because the numerical method constrains the wave 
function errors in a local fashion. 

Another group of test applications is comprised of well-studied leading order finite nuclear mass 
and relativistic corrections for the helium ground state. A straightforward computation reaches near 
state-of-the-art accuracy without requiring the implementation of any special-purpose numerics. 

All the relevant quantities tested in this paper - energy eigenvalues, S-state expectation values 
and bound-bound dipole transitions for S and P states - converge exponentially with increasing 
resolution and do so at roughly the same rate. Each individual calculation samples and weights the 
configuration space wave function uniquely but all behave in a qualitatively similar manner. Quan- 
tum mechanical matrix elements are directly and reliably calculable with pseudospectral methods. 

The technical discussion includes a prescription for choosing coordinates and subdomains to 
achieve exponential convergence when two-particle Coulomb singularities are present. The pre- 
scription does not account for the wave function's non-analytic behavior near the three-particle 
coalescence which should eventually hinder the rate of the convergence. Nonetheless the effect is 
small in the sense that ignoring the higher-order coalescence does not appear to affect adversely the 
accuracy of any of the quantities reported nor the rate at which errors diminish. 
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I. INTRODUCTION 



The aim of this work is to test and validate the pseu- 
dospectral method as a high-precision few-electron prob- 
lem solver, capable of calculating state-of-the-art preci- 
sion matrix elements. The helium atom has been studied 
extensively since the birth of quantum mechanics and so 
makes a great testbed problem. High-precision work con- 
tinues to this day to infer fundamental constants such as 
the fine structure constant (see Ref. |2|) and the electron- 
proton mass ratio (see Ref. 0]) by comparing theoret- 
ical and experimental measurements. Any theoretical 
method which may be applied to a variety of problems 
{e.g. high-precision relativistic corrections, different in- 
teraction potentials, excitation levels, symmetries, etc.) 
without tinkering with or modifying the basis and which 
has direct, rigorous control of local errors serves as a 
complementary approach to the variational method. 



Methods based on the variational principle, in which 
the expectation value of the Hamiltonian is minimized 
with respect to the parameters of a trial wave function, 
are the most widely used techniques for finding an ap- 
proximate representation of the ground state. The cal- 
culated energy is an upper bound to the exact energy 1 If 
one regards the best approximate wave function as first 
order accurate then the variationally determined energy 
eigenvalue is second order accurate. Small errors in the 
energy eigenvalue of a given state imply that the square 
of the wave function is accurate in the energy- weighted 
norm but it does not follow that local wave function er- 
rors are also small. In practical terms, while the vari- 
ational approach excels at determining energy eigenval- 
ues it does not generally achieve comparable accuracy in 
quantum mechanical matrix elements formed from the 
wave function. 

To achieve ever-more accurate energies and/or wave 
functions in the variational approach one must select a 
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1 The method is not limited to ground states. A trial wave func- 
tion, exactly orthonormal to all lower energy states, has calcu- 
lated energy which is an upper bound to the exact result for the 
excited state. 



sequence of trial functions capable of representing the 
exact solution ever-more closely. The choice of a good 
sequence entails more than a little art and intuition, es- 
pecially for a nonstandard problem where one may have 
only a vague idea what the ultimate limit looks like. A 
sequence of increasing basis size n may be said to con- 
verge exponentially if the errors are proportional to e~ an 
for some positive constant a. This most favorable out- 
come is achieved only if the basis can reproduce the an- 
alytic properties of the exact wave function. Otherwise, 
convergence is expected to be algebraic, i.e. oc n~ 2 , or 
worse. 

Recently, we applied pscudospectral methods to solve 
the nonrelativistic Schrodinger equation for helium and 
the negatively charged hydrogen ion with zero total angu- 
lar momentum [l|. We found exponentially fast conver- 
gence of most quantities of interest including the energy 
eigenvalues, local energy errors (e.g. (H^)/^> — E as a 
function of position) and Cauchy wave function differ- 
ences. Only the error in the logarithmic derivative near 
the triple coalescence point had discernibly slower conver- 
gence, presumably due to the logarithmic contributions 
located there 0-(y] • The key virtues of the pseudospectral 
approach were: no explicit assumptions had to be made 
about the asymptotic behavior of the wave function near 
cusps or at large distances, the Schrodinger equation was 
satisfied at all grid points, local errors decreased exponen- 
tially fast with increasing resolution, and no fine tuning 
was required. 

In this article, we extend our previous work to higher 
angular momentum calculations and utilize the results to 
evaluate matrix elements for combinations of states. To 
be systematic, we consider two sorts of matrix elements: 
the dipole absorption oscillator strength (between S and 
P states) and first-order mass polarization and a 2 rel- 
ativistic corrections to the nonrelativistic finite-nuclear- 
mass Hamiltonian (for the S ground state) . All have been 
the subject of extensive investigation. Our main focus is 
on testing the pseudospectral method's capabilities by re- 
calculating these quantities and comparing to effectively 
"exact" published results. 

The plan of the paper is as follows. The first four sec- 
tions are largely background: $11] provides an overview 
of the pseudospectral method; JjIIII describes the two- 
electron atom, the Bhatia-Temkin coordinate system, the 
expansion of the wave function in terms of eigenstates 
and the form of the Hamiltonian; £HVI defines length, ve- 
locity and acceleration forms for the oscillator strength 
and related sum rules. The next two sections detail our 
pscudospectral method of calculation and those read- 
ers primarily interested in seeing the results may skip 
to Will Wl gives a prescription for how to choose co- 
ordinates and subdomains for second order partial dif- 
ferential equations and outlines the special coordinate 
choices needed to deal with the Coulomb singularities. 
WII schematically describes how overlapping and touch- 
ing grids are coupled together and how symmetry is 
imposed on the wave function. Will presents the first 



group of test results on energies and oscillator strengths. 
The convergence rate of all quantities is studied in de- 
tail. Willi and i flXI review lowest-order corrections to 
the Hamiltonian due to finite nuclear mass and finite a. 
i jXl presents the second group of test results for individ- 
ual corrections to the ground state of He. £1X11 summa- 
rizes the capabilities and promise of the pseudospectral 
method. 

The appendix is divided into four parts. Appendix lAl 
gives the explicit form of the Hamiltonian operator used 
in this article. Appendix [B] describes how the Hamilto- 
nian matrix problem is solved, gives details of the eigen- 
value solver method, and how quantum mechanical ma- 
trix elements are calculated once the wave function is 
determined. Appendix [C] gives the particular equations 
for calculating the oscillator strengths and expectation 
values. Appendix [D] discusses and tabulates past work 
done to calculate oscillator strengths. 



II. REVIEW OF PSEUDOSPECTRAL 
METHODS 

Pseudospectral methods have proven success in solv- 
ing systems of partial differential equations germane to 
the physics in a wide variety of fields including fluid dy- 
namics Ml, g eneral relativity |8|, M\ > and quantum chem- 
istry [10l- |lS fl. Some problems in one-electron quantum 
mechanics [19L [20| have been treated but only recently 
has the method been applied to the case of fully corre- 
lated, multi-electron atoms [lj . Pseudospectral methods 
are discussed in some generality in Refs. [l|, l9l. l2l| - [24| . 

The pscudospectral method is a grid-based finite differ- 
ence method in which the order of the finite differencing 
is equal to the resolution of the grid in each direction. 
As the grid size increases it becomes more accurate than 
any fixed-order finite difference method. If a solution is 
smooth over an entire domain (or smooth in each sub- 
domain) the pscudospectral method converges exponen- 
tially fast to the solution. A spectral basis expansion 
and a pseudospectral expansion of the same order are 
nearly equivalent having differences that are exponen- 
tially small. 

The grid points in the pseudospectral method are lo- 
cated at the roots of Jacobi polynomials or their antin- 
odes plus endpoints. They are clustered more closely near 
the boundary of a domain than in its center. Such an 
arrangement is essential for the method to limit numer- 
ical oscillations sourced by singularities beyond the nu- 
merical domain [25| . These singularities typically occur 



in the analytic continuation of solutions to non-physical 
regimes and/or from the extension of coordinates beyond 
the patches on which they are defined to be smooth and 
diffcrcntiable. The grid point arrangement facilitates a 
convergent representation of a function and its derivative 
across the domain of interest. The interpolated function 
is more uniformly accurate than is possible using an equal 
number of equidistant points, as is typical for finite dif- 



ference methods. 

Consider the problem of the pscudospcctral represen- 
tation of an operator like the Hamiltonian. The full do- 
main is multi-dimensional but focus for the moment on 
a single dimension of the domain. Let {X. k }k=i.2,...N be 
the roots of an iVth order Jacobi polynomial enumerated 
by k. Let X stand for an arbitrary coordinate value in 
the dimension of interest. Define the one dimensional 
cardinal functions 



,v 



Q[x] = I 



fe=i 



x-x fc 
XJ - x fe 



and note the relation 



C 3 [X k 



(1) 



(2) 



follows. Now let the rid-dimcnsional grid be the tensor 
product of the individual, one dimensional coordinate 
grids labeled by X^ for i = 1 to n^. The corresponding 
cardinal functions are 



a property that may introduce non-physical effects, e.g. 
(X K \Hps\Cj) may possess complex eigenvalues. Gen- 
erally, unphysical artifacts quickly reveal themselves as 
resolution increases. An examination of the eigenvalue 
spectrum shows that the complex eigenvalues do not con- 
verge, permitting separation of physical and unphysical 
values. 



III. THE NONRELATIVISTIC 
TWO-ELECTRON ATOM 

Two-electron atoms are three-particle systems requir- 
ing nine spatial coordinates for a full description. In the 
absence of external forces, three coordinates are elimi- 
nated by taking out the center-of-mass motion. In the 
infinitc-nuclear-mass and nonrclativistic approximations 
the Hamiltonian is 



H = -\ip\+pl) + V, 



(7) 



c j[ x ] = Il^'wP^w]' 



(3) 



; = 1 



where subscript J = {j(i),j(2)i ■ ■ ■ ij(n d )} and unadorned 
X = {Xm), X( 2 ), ■ • ■ , X(„ d )}. These multi-dimensional 
Cardinal functions have the property 



Cj[X K ] = 6 



j • 



where the grid point X K = {X^ , X& , 



,X&}. They 



form a basis in the sense that a general function / can 
be written 



/[X] = £/[X J ]C./[X], 



(5) 



where /[X* 7 ] is a pseudospectral coefficient ("pseudo" be- 
cause it is more easily identified as the function value at 
the grid point). 

Let the position X K and cardinal Cj eigenstates be 
denoted \X K ) and \Cj), respectively. The pseudospectral 
approximation to the Hamiltonian is 



H PS =Y,\X K )(X K \H\Cj){C, l 



(0) 



JK 



where H is the full Hamiltonian operator. In practice, the 
matrix (X. K \Hps\Cj) is truncated and then diagonalizcd 
to find the energy eigenvalues. When the wave function 
is represented by a pseudospectral expansion the eigen- 
vectors are simply the function values at the grid points. 
In a spectral representation, by contrast, the eigenvec- 
tors are sums of basis functions. It is often more conve- 
nient and efficient to work with the local wave function 
values directly. On the other hand, the truncated op- 
erator Hps need not be Hcrmitian at finite resolution, 



where pi^ are the momenta of the two electrons and the 
potential is 



V 



z 



z 

I 2 



l 

ri2 ' 



(8) 



where Z is the nuclear charge, and r\, r%, and r\i are the 
magnitudes of the vectors pointing from the nucleus to 
each electron and of the vector pointing from one electron 
to the other, respectively Here and throughout this arti- 
cle, atomic units are used. For the infinite-nuclear-mass 
approximation, the electron mass is set to unity; for a 
finite nuclear mass, the reduced mass of the electron and 
nucleus is set to one. The fully correlated wave functions 
are six-dimensional at this stage. 

A further reduction is straightforward for S states. 
Hylleraas |26| proposed the ansatz that the wave func- 
tion be written in terms of three internal coordinates. 
Typical choices for these coordinates are r 1: r 2 , and Ty^. 
Alternatively, m may be replaced by #12, the angle be- 
tween the two electrons. The S state is independent of 
the remaining three coordinates that describe the orien- 
tation of the triangle with vertices at the two electrons 
and nucleus. 

The situation for states of general angular momentum 
is more complicated. Bhatia and Temkin |27| introduced 
a particular set of Euler angles {9, $, ^} to describe the 
triangle's orientation. They defined 2 a set of general- 
ized spherical harmonics D v tdm which are eigenstates of 
operators for the total angular momentum, its z compo- 
nent, total parity ({ri,r2J — > {— ri, — 1*2})) and exchange 



2 The symbols used here are slightly different than those of [27 
that the equations can be written in a simplified form. 
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^12-Dk( 



(-1) 
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K,lm' 



(9) 
(10) 

(11) 
(12) 



The superscript v takes on values v = and 1 while the 
integer subscript ft obeys < k < I. The quantum num- 
ber k is the absolute value of an angular momentum-like 
quantum number about the body-fixed axis of rotation. 
Even/odd ft determines the parity eigenvalue while the 
combination 1 + k+is determines the exchange eigenvalue. 
This basis is especially useful since each of the four op- 
erators above commutes with the atomic Hamiltonian, 
Hq . The spatial cigenfunction ipkims [ri , r 2 ] for total spin 
s, total angular momentum /, z-component of angular 
momentum m, and parity k = ±1 satisfies 



L 2 ip khns [ri,r 2 ] 

L z ipki m s[ri,r 2 ] 

n^«m S [ri,r 2 ] 

£i2^kims[ri,r 2 } 

Equations [9 1161 imply 



1(1 + l)ip klms [r 1 ,r 2 ] (13) 

rmpki ms [r 1 ,r 2 ] (14) 

kipklms[ri,r2] (15) 

(-l) s V^m S [ri,r 2 ]. (16) 



b 



klr 



1 <- 



(17) 
where the prime on the sum means that k is restricted 
to even (k = 1) or odd (fc = —1) numbers if parity is 
even or odd, respectively, and g v Kls is a real function of 
the internal coordinates. The convenience of the Bha- 
tia and Temkin [27| coordinate choice is most evident in 
how one imposes total antisymmetry of the wave func- 
tion. The spin singlet (triplet) must have a symmetric 
(antisymmetric) spatial wave function. The properties of 
the D^ lm functions reduce this requirement to 



£i29 v Kls = (-1) 



l/-\-K-\-l-\-S 



9kIs 



(18) 



The total antisymmetry of a wave function with given k, 
I, m and s follows by imposing the above requirement 
under r\ -s-> r 2 on each radial function for each v and 
ft. Note that (— l) K + i+s is fixed directly by the wave 
function's k, I and s. The same requirement applies to 
both singlet and triplet states up to the difference in the 
value of s. 

The full six-dimensional Schrodingcr equation for given 
I, s, even/odd parity, and any to yields I or I + 1 
(depending on these quantum numbers) coupled three- 
dimensional equations for g^. ls . The indices for g satisfy 
7 = or 1 and < ft < / with even or odd ft for even or 
odd parity. The equations are 



= (H S - E)gl ls 



l l 

EE 

i/=0n=-l 



HZ 



' n 

vtin-3 ti J r 2n,l,s ' 



(19) 



where Hs is the part of the Hamiltonian operator that 
survives for S states. The summation enumerates cou- 
plings with 7 7^ v and/or different k as well as terms 
that are intrinsic to non-S-states. 

Appendix [S] gives the explicit forms of the operators 
H s and iJJ K „. 



IV. REVIEW OF THE OSCILLATOR 

STRENGTH AND DIPOLE RADIATIVE 

TRANSITIONS 

The oscillator strength quantifies the coupling between 
two eigenstatcs of Hq on account of interactions with a 
perturbing electromagnetic field. It is fundamental for 
interpreting spectra, including the strength and width 
of atomic transitions and the lifetimes of atomic states. 
Sites generating spectra of interest are ubiquitous. They 
include earth-based laboratories, photospheres of the Sun 
and distant stars, and the near vacuum between the stars 
where traces of interstellar matter radiate. The specific 
applications of the oscillator strength are correspondingly 
diverse. For example, in laboratories the technique of 
laser spectroscopy is used to measure energy splittings 
and frequency-dependent photoabsorption cross sections 
of highly excited states. Knowledge of the transition 
probability matrices is needed to interpret which states 
have been directly and indirectly generated. The tran- 
sitions are driven by collisional and radiative processes, 
the latter given in terms of oscillator strengths. In an 
astrophysical context, on the other hand, observations 
of stellar emission require oscillator strengths for infer- 
ring chemical abundances from absorption or emission of 
radiation [28|, [29| • Oscillator strengths have widespread 
utility. 

The practical difficulty in calculating the oscillator 
strength value is the accurate representation of the initial 
and final wave functions. Almost from the very begin- 
ning of the development of quantum mechanics helium, 
having but two electrons, has served as a testing ground 
for new theoretical approaches. Appendix [D] presents a 
brief, schematic description of the rich history of such 
improvements in the service of oscillator strength calcu- 
lations. 

Following Baym [30( and Bethe and Salpeter [31|, the 
nonrelativistic Hamiltonian of a two-electron atom in 
the presence of an electromagnetic field (infinite-nuclear- 
mass approximation) is 



Hem — H n 



H 



int s 



(20) 



where Hq is the Hamiltonian for the isolated atom (Eq. 
[7]) and .Hint describes the interaction of the atom with 
radiation, 



-Hint = 






Pi ■ Aj + A,, ■ p, 
2c 



2c 2 



9: 



(21) 



where A^ and (pi are the vector and scalar potential, re- 
spectively, at the location of the ith electron (excluding 



the atomic Coulomb interactions included in V), and c is 
the speed of light. If the photon number density is small 
then the second term, corresponding to two-photon pro- 
cesses, is much smaller than the first and if one adopts 
the transverse gauge then the third term is zero. With 
these assumptions the non-zero terms are the ones linear 
in the vector potential. 

Only electric dipolc-mcdiatcd transitions and the as- 
sociated /'s are considered in this article. The length, 
velocity and acceleration forms for the oscillator strength 
[Si are 



respect to the ground state. The rules [3J, |35| include 



■' i j 
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f" 



-(Ej- £,-)|<jTOI 2 
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\(J I A|i> 



(22) 
(23) 
(24) 



Here Ei and Ej are the energies of the initial and final 
states. The two-particle operators are 



R = n + r 2 


(25) 


P = Pi +P2 


(26) 


A - ZTl Zr2 


(27) 



i.e. the position, momentum and acceleration electron 
operators. Appendix [C] presents explicit expressions for 
/ used in the calculations. 

If the wave functions, energies, and operators were ex- 
act, all three forms would give identical results. How- 
ever, in a numerical calculation the agreement may be 
destroyed whenever the operator commutator rule 



P = i[H Q ,R} 



(28) 



is violated. Approximations to the operators (Ho, P, 
or R) and to the initial and final eigenstates are possi- 
ble sources of error. Good agreement between the three 
forms at a fixed resolution has sometimes been taken to 
be an indication of an accurate answer. Such agreement 
is ultimately necessary as resolution improves but the 
closeness of the agreement is insufficient to infer the ac- 
curacy at a fixed resolution J32|, [33| . A more stringent 
approach involves two steps: first, for each form check 
that the matrix element converges with resolution or ba- 
sis size and, second, that the converged answers for dif- 
ferent forms agree. 

The oscillator strengths fo n for transitions, 1 1 S — > n l Y 
of helium obey a family of sum rules. For integer k define 



S(k) = J2\ AE °"\ k f0r, 



S(-l) = 


3<(n+r 2 ) 2 ), 


5(0) = 


2, 


5(1) = 


4 , - . 

-g(#0 "Pi -P2), 


5(2) = 


M(o( ri ) + <5(r 2 )), 



(30) 

(31) 

(32) 
(33) 



where the expectation values on the right hand side refer 
to the ground state. 

In principle, these sum rules provide consistency checks 
on theoretically calculated oscillator strengths. However, 
the explicit evaluation of S(k) (Eq. [23) is difficult. Mul- 
tiple methods are needed to handle all the final states, 
which include a finite number of low energy highly corre- 
lated states, a countably infinite number of highly excited 
states, and an uncountably infinite number of continuum 
states. Ref. [36( inferred that the two sides of Eqs. I3"01l3"3l 
agree to about one percent based on a combination of the 
most reliable theoretical and/or experimental values for 
fon- 

This article exemplifies the capabilities of the pseu- 
dospectral approach by evaluating the l-'-S — > i^Y oscil- 
lator strength, a physical regime in which strong electron 
correlations are paramount, and a set of expectation val- 
ues for operator forms, some of which appear on the right 
hand side of the sum rules. 



V. VARIABLES AND DOMAINS 

This section details an important element of the appli- 
cation of the pscudospcctral method: the choice of coor- 
dinates and computational domains. 

To achieve exponentially fast convergence with a pseu- 
dospectral method, it is imperative that the solution be 
smooth. The presence of a singular point may require a 
special coordinate choice in the vicinity of the singularity 
or a different choice of effective basis. Handling multiple 
singularities typically requires several individual subdo- 
mains, each accommodating an individual singularity. It 
is useful to have a guide for choosing appropriate coordi- 
nates. 

The ordinary differential equation 



(— 

\dX 2 



Pa [X] d q a [X] 



x-adx {x- a y 



/ = o 



(34) 



(29) 



with p a [X] and q a [X] analytic at X = a has a regular 
singular point at X = a. The basic theory of ordinary 
differential equations (ODE's) [331 states that / has at 
least one Frobenius-type solution about X = a of the 
form 



where the summation is over all P states, including the 
continuum. Here, A.Eq„ is the energy difference with 



f[X] = (X-a) t «J2c n (X-a) n , 



(35) 



71 = 



where the coefficients c n can be derived by directly plug- 
ging into Eq. [M] and t a is the larger of the two solutions 
to the indicial equation 



ta{ta - 1) +Pa[a]t a + q a [a] = 0. 



(36) 



Exponential convergence of the pseudospectral method 
for a differential equation of the form of Eq. |3U requires 
t a be a non-negative integer. This must hold at each 
singularity a in the domain (as well as all other points). 3 
A simple example is the Schrodinger equation for 
a hydrogenic atom expressed in spherical coordinates 
{Xi, X2, X3} = {r, 9, </>}. The radial part of the full wave 
function -R„z[r] satisfies 



dr 2 



2 d 
r dr 



1(1 + 1) - 2Zr - 2Er 2 



Rni = 0- (37) 



A comparison with Eq. [34] yields po [0] = 2 and qo [0] = 
—1(1 + 1), which gives to = I, the well known result for 
hydrogenic wave functions. The reduction of the partial 
differential equation (PDE) into an ODE having non- 
negative integer to tells us that spherical coordinates are 
a good choice for solving hydrogenic wave functions using 
pseudospectral methods. A bad choice would be Carte- 
sian coordinates {Xi,X2,X3} = {x, y, z}. The ground 
state has the form 



i> oc e ~ z V* 2 +y 2 +z 2 



(38) 



This solution has a discontinuity in its first derivatives at 

x = y = z = 0: 

lim — ^ lim — . (39) 

x,y,z^0+ OX,y, Z x,y,z-*0~ ox, y, z 

Other solutions have a discontinuity of first or higher 
derivatives at the same point. The pseudospectral 
method would not handle these well and convergence 
would be limited to being algebraic. 

An arbitrary second order PDE may have singulari- 
ties that occur on complicated hypersurfaces of different 
dimensionality. Deriving the analytic properties of a so- 
lution near such a surface is a daunting task. The general 
idea is to seek a coordinate system such that the limit- 
ing form of the PDE near the singularity looks like an 
ODE of the sort that pseudospectral methods are known 
to handle well. 

For example, in a three-dimensional space, assume the 
singularity lies on a two-dimensional surface. First, seek 
a coordinate system such that the surface occurs at Xi = 



3 The full class of one dimensional problems for which pseudospec- 
tral methods converge exponentially fast is larger than this de- 
scription. The method needs the solution to be smooth which is 
a weaker statement than that it be analytic. This distinction is 
not material for the singular points discussed here. 



a. 4 Second, focusing on Xi, seek coordinates so that is 
possible to rewrite the PDE in the form 



ifi 



Pa[X] 8 



Qa\X] 



dXi 2 Xi-oaxx (Xi-o) 



/ = (40) 



where P a and Q a are linear second order differential op- 
erators that do not include derivatives with respect to 
Xi. Finally, seek coordinates such that P a and Q a are 
analytic with respect to Xi at a. 

Unfortunately, even if one succeeds in finding such a 
coordinate system, the theorem of ODEs docs not gener- 
alize to PDEs, i.e. there is no guarantee that / is analytic 
near a. A celebrated example is exactly the problem 
of concern here, i.e. the Schrodinger equation for two- 
electron atoms. Three coordinates are needed to describe 
the S state. In hyperspherical coordinates ({Xi, . . . } = 
{p, . . . } where p = \/t\ + r 2 ), Schrodinger's equation 
matches the form of Eq. [40] for Xi = p and a = 0. 
This is the triple coalescence point, a point singularity in 
the three-dimensional subspacc spanned by the coordi- 
nates n , r<i , and r\% . The electron- nucleus and electron- 
electron singularities (two-body coalescence points) are 
one-dimensional lines in this subspacc that meet at p = 0. 
Bartlett |4| proved that no wave function of the form 



v> = f>r 



(41) 



?i=0 



where A n is an analytic function of the remaining vari- 
ables will satisfy the PDE. Fock's form for the solution 
ills 



oo L«/2J 

v^ = E E B nmP n (\og P y 



(42) 



n— 771—0 



where B nm is an analytic function of the remaining vari- 
ables. The presence of the logp terms in the wave func- 
tion is an important qualitative distinction between a 
solution having two- and three-body coalescence points. 
Some properties of the solution near p = have been 
reviewed in our previous article [l|. For example, Myers 
et al. [38} showed that the logarithmic terms allow the lo- 
cal energy (Hip) /ip near p = to be continuous. Despite 
this property, they have only a slig ht effect on the con- 
vergence of variational energies [391 ] . By many measures 
of error the triple coalescence point does not affect pseu- 
dospectral calculations until very high resolutions [l| . 



4 A zero- or one-dimensional singularity can be made to look two- 
dimensional by a coordinate transformation. For example, in 
the previous example, which has used spherical coordinates, the 
Coulomb singularity appears at r = 0. This point is approached 
on a two-dimensional sphere of constant radius by taking the 
limit as a single coordinate, the radius, approaches zero. 



As a point of principle, however, no simple coordinate 
choice can hide the problems that occur at the triple 
coalescence point, and no special method for handling 
this singularity is given here. Elsewhere (p ^ 0) our rule 
of thumb is the following: coordinates are selected so that 
the singularity may be described by X% — a with P a and 
Q a satisfying 



P a = £(X,-a) B j 



n=0 



Q a = ^2(Xi-a) n q an , 

n=0 



(43) 
(44) 



in a neighborhood about Xj = a. Here, p an and q an 
are linear differential operators not containing X, or its 
derivatives. 

The singularities of the Hamiltonian, given in detail in 
Appendix [2 are of two types. The physical singularities 
at n, r 2 , and rvi = were explored in Ref. |l|. One of 
the essential virtues of hypersphcrical coordinates is that 
p 7^ implies these coalescences have separate neighbor- 
hoods. Therefore, the prescription is to seek separate 
coordinates satisfying eqs. 03] and |U] in the vicinity of 
each singularity. 

There are also coordinate singularities at 6\ 2 = and 
7r which correspond to collincar arrangements of the two 
electrons and nucleus. These singularities were com- 
pletely absent in our previous treatment of S states |l| 
where C — — cos #12 and B = — cos ji\ 2 (fin is defined be- 
low) were the third coordinates in different subdomains. 
Now, to accommodate the singularities' presence in the 
Hamiltonian for general angular momentum make the 
slight change to use #12 and j3\ 2 instead. 

Starting with the internal coordinates n, r 2 and #12 
one defines p, (f>, £, and x by 
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FIG. 1: (Color online). This is the arrangement of grid points 
of the three domains at a constant value of p in and #12 co- 
ordinates for n — 20. Note that the point density becomes 
larger at the boundary of each subdomain and that no grid 
points sit on the Coulomb singularities. The blue circles, red 
crosses, and green pluses belong to domains D\, D2, and D3, 
respectively. D\ and D2 are rectangular domains, while D3 
has the curved boundary in <j>, #12 coordinates but is rectan- 
gular in f, /3i2 coordinates. The electron-proton singularity 
occurs on the left side (solid line at cj> = 0). The entire line 
corresponds to one physical point. The electron-electron sin- 
gularity occurs at the lower right hand corner (solid disk at 
<t> — 7r/4, 812 = 0). A line of symmetry falls on the right side 
(dashed line at <j> — n/4 where n = r^)- 



Eqs. [43j and |44] arc satisfied by selecting 

{Xi,X 2 ,X 3 } = {x,(j>,6i 2 } or {x, C/812} in three sepa- 
rate domains 



''1 


= p cos 4> 




(45) 


r-2 


= p sin <p 




(46) 


'12 


= p^/2si^\( > 




(47) 


V^sinC 


= y/l — cos9\ 2 sin 2<p 


(48) 


COS /3i2 


cos 2(p 




(49) 


V 1 — cos 2 9\2 sin 2 


20 


X 


1-p 

l + p 




(50) 


The full ranges of these variables are 






\n 


< ri,r 2 ,p < 00 

-r 2 \ < r 12 < n + r 2 

< 012, 012 < 7T 

0< 0,C<tt/ 2 

-1 < X < 1. 




(51) 



The purpose of coordinate x is to map the semi-infinite 
range of p to a finite interval. 



A : 


-1< x < 1, 


D 2 : 


-l<x< 1, 


D 3 : 


-1<X< 1, 



0< <j>< 



< 



i_ 

2' 

< £ 

— it' 



o<C< 



-1 < COS 6<i2 < 1 
-1 < C0S6»12 < § 
-1 < COS ^12 < 0, 



(52) 

spanning only half the space defined by the inequali- 
ties (|51|) due to the symmetry in the Hamiltonian about 
7"i = r 2 . Fig. Q] illustrates the layout of the three do- 
mains at fixed p. The coordinate systems in domains D\ 
and D3 were developed to handle the electron-proton and 
electron-electron singularities, respectively. The choice 
of coordinates in domain D 2 was more arbitrary, and for 
simplicity was chosen to be the same as in domain D±. 
This particular choice allows for no overlap between do- 
mains D\ and D 2 and makes the symmetry condition 
(Eq. [T8"|) at n = r 2 , <f> — 7r/4, or fii 2 = ir/2 easy to ap- 
ply. The remaining electron- nucleus singularity, n = 0, 
is implicitly accommodated by the spatial symmetry of 
the wave function. The three domains must jointly de- 
scribe the full rectangle but the specific choice for edges 
at0 = £ = l/2is arbitrary. 



VI. BOUNDARY CONDITIONS 
A. Internal boundary conditions 

It is necessary to ensure continuity of the wave function 
and its normal derivative at internal boundaries. There 
are two ways in which the subdomains can touch: they 
can overlap or they can barely touch. For clarity, con- 
sider a one-dimensional problem with two domains. Let 
the first domain be domain 1 and the second be domain 
2 with cxtrcma Xi iinin < X 2 , m in < X 1>max < X 2 , max , 
where the 1 and 2 refer to domain number. The first 
case corresponds to X2, m i n < Xi max and the second to 
X2,min = Xi >max = X» . For both cases, exactly two con- 
ditions are needed to make the wave function and its 
derivative continuous. The simplest choice for the first 
case is 



■*M X l,max] = ^2[Xl,I 
01 [X 2 ,min] = V^P^r 

and for the second case is 



0i [x„] = MM 



dx 



?/>i[X* 



dX 



M**\- 



(53) 
(54) 



(55) 
(56) 



For multi-dimensional grids, the situation is analogous. 
The conditions are applied on surfaces of overlap. In 
this case the derivatives are surface normal derivatives or 
any derivative not parallel to the boundary surface. On 
a discrete grid, a finite number of conditions are given 
which, in the limit of an infinitely fine mesh, would cover 
the entire surface. Additional discussion and illustrations 
of the technique are in Ref. [l[ . 



B. The symmetry condition 

The Hamiltonian (see appendix [A| is symmetric with 
respect to particle exchange {r\ <-> r%). Therefore, there 
are two types of eigenstates: those with symmetric spa- 
tial wave functions (singlets) and those with antisym- 
metric spatial wave functions (triplets). The radial wave 
functions g v Kls satisfying the appropriate symmetry must 
obey Eq. [18j More explicitly 



d9li s 



9nU 



<>u',: 



4=7174 



9/3i 



/3l2 = 7T/2 



if £ is even 



j. a = <7^i.l o m if £ is odd 



(57) 



where £ = v + K + l + s. 



VII. ENERGY AND OSCILLATOR STRENGTH 
RESULTS 

This article generalizes the pseudospectral methods 
previously developed for S states to the general angular 
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FIG. 2: (Color online). The logarithm base 10 of the energy 
error (AQ) of both the lowest energy S state and P state of 
helium. The dark blue circles are for the 1 X S state and the 
light red crosses for the 2 J P state with dashed blue and dotted 
red fits, respectively (see Tab. [I]). 



momentum case, calculates oscillator strengths for tran- 
sitions, and tests how different measures of wave function 
errors vary with resolution. 

The most widely quoted number to ascertain conver- 
gence is the energy which gives a global measure of ac- 
curacy. Figure [5] shows the energy errors for the 1 X S 
and 2 1 P states of helium. Here and throughout the re- 
sults sections the high precision values of Drake [35( are 
taken to be exact. The energy error for both states de- 
creases exponentially with resolution. Convergence for 
the S state is similar to that reported in Ref. [l[ with 
slight differences related to a different choice of coor- 
dinates. The current calculation extends to basis size 
n = 23 for S states and n = 20 for P states instead of 
n = 14 for only S states in Ref. [1(. 

A common feature of the energy convergence and all 
other convergence plots in this article is non-monotonic 
convergence. This method is not variational, so there is 
no reason to expect monotonic convergence. Calculated 
quantities can fall above or below their actual value, with 
error quasi-randomly determined by the exact grid point 
locations. The jumps decrease in magnitude as the reso- 
lution is increased. 

As described in Sec. IIV1 there are three commonly 
used forms for the oscillator strength. The length, veloc- 
ity, and acceleration forms depend most strongly on the 
value of the wave function at positions in configuration 
space corresponding to large, medium, and small sepa- 
rations. Sometimes the relative errors are used to infer 
where the wave function is more or less accurate. It has 
been observed that for most variational calculations, the 
acceleration form tends to be much less accurate than the 
other two forms, suggesting errors in the wave function 
at small separation that have little effect on the varia- 
tional energy. The length and velocity forms give results 
of roughly comparable accuracy. 

The oscillator strength of the 1 X S — > 2 1 P transition 
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FIG. 3: (Color online). The logarithm base 10 of the error 
(AQ) in the oscillator strength of the 1 1 S — > 2 X P transition 
of helium. The dark blue circles are for the length form, the 
light red crosses are for the velocity form, and the green pluses 
for the acceleration form with dashed blue, dotted red, and 
dot-dashed green fits, respectively (see Tab. [I]). 



FIG. 4: (Color online). The logarithm base 10 of the error 
(AQ) in the average of the length, velocity, and accelera- 
tion forms of the oscillator strength (dark blue circles) and 
their standard deviation (light red crosses) for the IS — > 2*P 
transition of helium, with dashed blue and dotted red fits, 
respectively (see Tab. [I]). 



was calculated using all three forms and Fig. [3] displays 
the errors. Here, all three forms give roughly the same 
results. At most resolutions the points lie nearly on top 
of one another and their fits are indistinguishable, indi- 
cating the wave function errors for small, medium, and 
large separations have roughly equal contributions to the 
numerically calculated oscillator strength. This may be 
due to the pseudospectral method's equal treatment of 
all parts of configuration space. 

It should be noted that the value used as the exact 
value [35| is given to seven decimal places. Consequently, 
the errors inferred for the highest resolution calculations 
in Fig. [3] are not too precise. There is little practical 
need for additional digits since a host of other effects in- 
cluding finite nuclear mass, relativistic, and quadrupolc 
corrections would confound any hypothetical, experimen- 
tal measurement of the oscillator strength to such high 
precision even if a perfect measurement could be made. 
Actual experiments struggle to obtain two percent preci- 
sion [4(|, an error larger than these effects. 

As pointed out by Schiff et al. [33| and reviewed by 
Hibbert [32(, the assumption that using the differences 
between the oscillator strength values from the differ- 
ent forms as a measure of the accuracy is not valid. 
Agreement is necessary but not sufficient. They suggest 
comparing calculated and extrapolated values. This lat- 
ter procedure is not straightforward for a pseudospectral 
method with non-monotonic convergence. We present a 
similar suitable check. Fig. |4]shows the average and stan- 
dard deviation of the error for the three forms as a func- 
tion of resolution. The standard deviation is about an or- 
der of magnitude (with a large scatter about that factor 
of ten) less than the average error at low and moderate 
resolutions but the trend lines suggest that the standard 
deviation may be approaching the average at the higher 
resolutions. A possible explanation is that the calcula- 



TABLE I: The fit parameters to all the convergence plots of 
quantities Q in this section. 



Q 


Figure 


A 


P 


£(1 X S) 


m 


2.5 x 10" 9 


0.40 


£(2 1 P) 


m 


5.2 x 10~ 9 


0.42 


J12 


El 


8.4 x 10~ 8 


0.40 


rv 
J12 


m 


9.2 x 10~ 8 


0.39 


fl2 


m 


8.6 x 10~ 8 


0.40 


ravg 
J 12 


® 


8.7 x 10" 8 


0.40 


fSD 

J12 


m 


2.2 x 10~ 8 


0.34 



tion at the highest resolutions is starting to become sensi- 
tive to the wave function truncation (see appendix IB 21) . 
This destroys the expected equality between the forms 
and each form converges to its own incorrect asymptotic 
value. The individual errors and the standard deviation 
become comparable. So at n = 20, we assume the stan- 
dard deviation and total error are equal and get a value 
for the oscillator strength of 0.27616499(27) which com- 
pares favorably to Drake's 0.2761647 [351 ] - 

All convergence data were fit to functions of the form 
AQ = A x 10~ /3 (™ _2 °) using the same procedure as in 
Ref. [lj. Because of uncertainty in the errors for the 
largest resolutions (n = 19 and n = 20) these points 
were not used in the fits of f[ 2 , /" 2 , fi2: an d fi^ 9 - The 
j3 parameter, which corresponds to the slope of the fits 
in the convergence graphs is roughly the same for all fits, 
with the exception of the standard deviation of the os- 
cillator strength forms. This behavior is consistent with 
our discussion of errors in the previous paragraph. 
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VIII. CORRECTIONS TO THE HAMILTONIAN 

Two small parameters appear in the full physical 
Hamiltonian: the ratio of the reduced mass of the 
electron- nucleus pair to the nuclear mass, n/M = 
1.37074563559(58) x 10" 4 jHIH (for 4 He) and the fine 
structure constant a = 7.2973525376(50) x 10" 3 [Mill- 
Here, the lowest order corrections in n/M and a are con- 
sidered. For very high-precision work, one needs the pcr- 
turbative corrections in powers of each small quantity. 

A. Finite nuclear mass correction 

The nonrelativistic (a ) Hamiltonian for two-electron 
atoms is 



H nl — H + H c 



H„ 



(58) 



where Ho is the fixed-nucleus approximation to the 
Hamiltonian with the electron mass set to /i, H cm is the 
kinetic energy of the center of mass, and H mp is the mass 
polarization term: 



H () 



1 



(pI+pD + v 



lie 



H„ 



2(M- 
7^ Pl 



2m e ) 



rem 



P2, 



(59) 
(60) 
(61) 



where V is the potential energy operator, m e is the elec- 
tron mass, p cm is the momentum operator of the center 
of mass, and reduced mass atomic units (/i = 1) are be- 
ing used. The second term is removed in center-of-mass 
coordinates and the last term provides the dominant non- 
trivial correction for finite nuclear mass (the trivial one 
being the scaling of the energy by m e / jjl). 



The separate contributions to the Hamiltonian are the 
mass- velocity (mass), two-body Darwin (D), spin-spin 
contact (SSC), orbit-orbit (00), spin-orbit (SO), spin- 
other-orbit (SOO), and the spin-spin (SS) terms. These 
are explicitly given by 



-^ mass 

H D = 

Hssc = 

Hoo = 

H so = 



«V 4 

% 

87ra 2 



><j 



(si • s 2 )<5(ri 2 ) 

a 2 /^Pi -P2 , ri 2 (ri2 • pi) • p 2 
2 



a 2 Z 



E 



''12 

T • s 



'12 



Hi 



SOO 



= 4E ^XPi •(8i + 2s i ) 



. r 



Hss — 



^( s i- s 2-^( s i' r i2)(s2-r 12 ; 

'12 \ '12 



(65) 
(66) 

(67) 
(68) 

(69) 

(70) 
,(71) 



where i and j can be 1 or 2, p^ and r^ are the momentum 
and position of the ith electron with respect to the nu- 
cleus, respectively, ri 2 is the vector pointing from the first 
electron to the second, and §i and 1^ are the one-electron 
spin and angular momentum operators of the ith elec- 
tron, respectively. The last three Hamiltonian terms are 
zero for 1 S states due to symmetry considerations. 

There are many higher order terms (see Refs. [2, H - 
|46| ) but these are not considered here. 



IX. MASS POLARIZATION AND 
RELATIVISTIC CORRECTION CALCULATIONS 



B. Relativistic corrections 

The Schrodingcr equation is a nonrelativistic approxi- 
mation to the true equation of motion. The lowest order 
relativistic corrections enter at order (a 2 ), as summarized 
in Ref. [43[ and repeated here. Note, all references in this 
article to orders in a are in Rydbergs. The Brcit-Pauli 
Hamiltonian encapsulates the correction 



Hbp = H n 



H IC 



(62) 



where -ff nr is the usual nonrelativistic Hamiltonian used 
in Schrodingcr 's equation and H TC \ is the lowest order 
relativistic correction. The latter can be further di- 
vided into non-finc-structure (NFS) and fine-structure 
(FS) contributions: 

#NFS = -ffmass + #D + -^SSC + H (63) 

#fs = Hso + Hsoo + Hss- (64) 



The mass polarization and low order relativistic correc- 
tions to the nonrelativistic Hamiltonian have been known 
for some time |3l|. The main challenge in calculating 
these terms is finding adequate unperturbed wave func- 
tions. Early calculations p7H50| were critical for com- 
paring experimental and theoretical energies, confirming 
that Schrodingcr's equation is correct in the nonrelativis- 
tic limit for helium. 

The develo pme nt of computers enabled Pekcris and 
coworkers |5lTi53l | and others [54M58J to reach theoretical 
uncertainties in the energy of about 10 -2 cm -1 . Such 
precision and the resulting precision in the wave func- 
tion allowed Lewis and Serafino [57| to calculate the fine 
structure constant from experimental measurements of 
the 2 3 P splitting. They obtained a" 1 = 137.03608(13) 
with an estimated uncertainty only surpassed at the time 
by the measurements of the electron anomalous magnetic 
moment (g — 2) (by a factor of two) and the ac Joscphson 
experiments (by a factor of four). 
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Drake and collaborators [4J, |59U64| and Pachucki and 
collaborators j2|, I65l - l74 | have pushed relativistic correc- 
tions for regular helium up to order a 5 and beyond us- 
ing a Hylleraas [26| type basis. Drake [63[ matched 
theoretical and observed energy differences in the J = 
0, 1 splitting of the 2 3 P state and determined a -1 = 
137.0359893(23). Drake cited a difference with the g - 2 
result 137.0359996(8) but agreement with the ac Joseph- 
son result 137.0359872(43) [63( . However, a similar calcu- 
lation of his using the observed J = 1, 2 splitting gives an 
unreasonable value [63J . Pachucki and collaborators have 
resolved the issue by finding errors in a 5 terms and by in- 
creasing the error estimate due to a 6 terms. Their most 
recent determination is or 1 = 137.03599955(64)(4)(368), 
where the first error is experimental, the second numer- 
ical, and the third is their estimated error from higher 
order terms [2| • This value agrees with the latest g — 2 
results but is not as precise [2|. 

An alternative approach is to use an even simpler basis, 
with surprisingly accurate results. Korobov and collabo- 
rators have used an exponential basis (see Refs. j75l[76|) 
to calculate very precise heliu m I77H82J (up to order a 4 ) 
and anti-protonic helium [3. I83l - l89j (up to order a 5 ) elec- 
tronic energies. The latter calculations have been used 
for the CODATA06 [U 53 recommended value of the 
electron-to-(anti)proton mass ratio. 



X. EXPECTATION VALUES 
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FIG. 5: (Color online). The logarithm base 10 of the error 
(AQ) in the expectation values of operators that scale as p 2 
for helium. The dark blue circles are for (r 2 ), the light red 
crosses are for (rf 2 ), and the green pluses are for (nr2 cos 612) 
with dashed blue, dotted red, and dot-dashed green fits, re- 
spectively (see Tab. 111)) . 



The aim of this section is to test the pseudospectral 
method's ability to represent the wave function in dif- 
ferent parts of configuration space and to compare the 
convergence rates of the errors with that of the ener- 
gies and oscillator strengths. For a representative set of 
calculations consider the expectation values of the oper- 



ators needed for leading order relativistic (Sec. IVIII B() 
and finite nuclear mass (Sec. IVIII A[) corrections, for the 
oscillator strength sum rules (Eqs. I3TJ1I55]) . interparticle 
distances, (V), and (V 2 ). These expectation values test 
different parts of the wave function as well as different 
types of operators. They are organized by the weighting 
of the wave function and used to draw inferences about 
local errors. 
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FIG. 6: (Color online). The logarithm base 10 of the error 
(AQ) in the expectation values of operators that scale as p 
for helium. The dark blue circles are for (ri) and the light 
red crosses are for (ru) with dashed blue and dotted red fits, 
respectively (see Tab. 111)) . 

Figure [5] displays results for expectation values related 
to sum rule S(—l) (Eq. [30]), i.e. quantities scaling like p 2 . 
These calculations are somewhat more sensitive to the 
wave function at large separation than, say, the normal- 
ization integral. In addition, they focus on parts of coor- 
dinate space which have low resolution compared to the 
coverage near the singularities. High accuracy is found 
for all three cases. 

Figure [5] displays results for expectation values of op- 
erators scaling like p similar to the length form of the os- 
cillator strength. Higher accuracy is obtained here than 
for the oscillator strength at equivalent resolutions. This 
can be explained by the smaller length scale set by the 
higher energy of the P state, which enters only into the 
oscillator strength calculations. So a greater resolution 
is needed for the same accuracy. 

Figure [7] displays results for expectation values related 
to the potential energy of charged particles, i.e. quan- 
tities scaling like 1/p. This probes the treatment of the 
singularities. The high degree of accuracy is evidence 
that these singularities have been treated correctly. 

Figure [8] displays results for expectation values related 
to the square of the potential energy, i.e. quantities scal- 
ing like 1/p 2 . These operators emphasize the singularities 
even further. One may expect that at a high enough in- 
verse power of p that the effect of the Fock logarithm 
become important and slow down convergence, but no 
evidence of that effect is apparent. 

Even the expectation values of delta functions, related 
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FIG. 7: (Color online). The logarithm base 10 of the error 
(AQ) in the expectation values of operators that scale as 1/p 
for helium. The dark blue circles are for (1/Vi) and the light 
red crosses are for (I/7-12} with dashed blue and dotted red 
fits, respectively (see Tab. Hi)) . 



FIG. 9: (Color online). The logarithm base 10 of the error 
( AQ) in the expectation values of delta function operators for 
helium. The dark blue circles are for (<5(ri)) and the light red 
crosses are for (<5(ri2)) with dashed blue and dotted red fits, 
respectively (see Tab. HT)) . 




FIG. 8: (Color online). The logarithm base 10 of the error 
( AQ) in the expectation values of operators that scale as 1/p 2 
for helium. The dark blue circles are for (1/ri), the light red 
crosses are for (l/r 2 2 ), the green pluses are for (l/rir^), and 
the black stars are for {l/rir-12} with dashed blue, dotted red, 
dot-dashed green, and solid black fits, respectively (see Tab. 

E}. 



to form the appropriate operators) appear to be just as 
accurate as the function values, even when they are most 
strongly weighted close to the electron-electron cusp, as 
is the case for the orbit-orbit interaction. 

The exponential rate of convergence and the magni- 
tude of the errors are roughly the same in all the calcula- 
tions of expectation values in Figs. IfjIfTUl This is reflected 
in the fits (sec Tab. [TTJ) . 
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to sum rule S(2) (Eq. [22]), the Darwin term H D (Eq. I55|) . 
and the spin-spin contact term H^sc (Eq. I67[). which are 
most sensitive to the Kato cusp conditions |9(| have the 
same convergence properties (See Fig. ^. This provides 
evidence that our choices of coordinates allowed the pseu- 
dospectral method to deduce and represent the solution 
in the vicinity of a cusp. It also shows that if one can han- 
dle the non-analyticities of the matrix element by hand, 
as is possible for delta functions (see appendix IB 3|) . one 
can still have exponentially fast convergence. 

The error in the mass polarization if mp (Eq. I61[) . used 
for the finite-nuclear mass correction and the calculation 
of the sum rule 5(1) (Eq. I52"j) . and the orbit-orbit terms 
Hoo (Eq. I68[) . i.e. quadratic momentum contributions, 
are shown in Fig. [TO] Calculations of derivatives (needed 



FIG. 10: (Color online). The logarithm base 10 of the error 
(AQ) in the expectation values of the mass polarization and 
the orbit-orbit interaction operators for helium. The dark 
blue circles are for (pi • P2) and the light red crosses are for 
{Hoo)/a 2 with dashed blue and dotted red fits, respectively 
(see Tab. [II]). 

These errors decrease until they reach roughly the level 
of error produced by truncating the wave function (see 
Sec. [B]) at the highest resolutions. The only easily dis- 
cernible differences are at low resolution for which the 
representation of the wave function at large p is certainly 
poor. It is unsurprising that the expectation values that 
scale as p 2 and p have larger errors at low resolution due 
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TABLE II: The fit parameters to all the convergence plots of 
quantities Q in this section. 



Q 


Figure 


A 


P 


(r?> 


m 


3.3 x 10" 11 


0.54 


(rh) 


m 


1.2 x 10~ 10 


0.53 


(ri -r 2 ) 


m 


4.9 x KT 11 


0.47 


(ri) 


m 


1.2 x 10" 11 


0.48 


<ri2> 


El 


7.5 x 1CT 11 


0.46 


(l/n) 


m 


1.1 x 1(T 10 


0.37 


(l/ria) 


m 


1.2 x 10" 10 


0.37 


<l/r?> 


m 


7.1 x 1(T 10 


0.36 


(l/r-? 2 > 


m 


3.1 x 10" 10 


0.38 


(1/rira) 


m 


3.9 x 10" 10 


0.37 


(1/nria) 


m 


2.6 x 1(T 10 


0.37 


(S(ri)) 


m 


2.4 x 10" 10 


0.36 


(S(ri2)) 


El 


6.5 x 10" 11 


0.36 


(Pi ■ P2> 


M 


2.4 x 10~ 10 


0.39 


(Hoo)/a 2 


M 


4.3 x 1(T 10 


0.38 



to the scarcity of points in the asymptotic tail of the wave 
function. 

All convergence data were fit to functions of the form 
A x 10 _,3 (" -23 ) using the same procedure as in Ref. [l|. 
The fit parameters are shown in Tab. fll] The most strik- 
ing feature is how similar the magnitudes of the errors are 
at n = 23. Also, the exponential parameter j3 is roughly 
the same for all expectation values and the energies and 
oscillator strengths (see Tab. HJ with the differences al- 
ready discussed. Indeed, as one increases resolution one 
increases the accuracy of all expectation values or oscil- 
lator strengths by roughly the same amount. 

The contributions to the total energy of the ground 
state of 4 He are summarized in Tab. IIIII The values 
from both this work and Drake's [35[ are given. For a 
wave function with a much lower precision in its eigen- 
value (nine decimal places compared to fifteen), nearly 
the same precision is obtained for the corrections to this 
eigenvalue. 



XI. CONCLUSIONS 

We developed a general prescription for choosing coor- 
dinates and subdomains for a pseudospectral treatment 
of partial differential equations in the presence of phys- 
ical and coordinate-related singularities. This prescrip- 
tion was applied to Schrddingcr's equation for helium to 
determine the fully correlated wave function. The treat- 
ment accounts for two-body but not three-body coales- 
cences. Other problems with Coulomb singularities can 
now be tackled with this method. 

We explored the fidelity of the pseudospectral 
method's results. The method attained exponentially 



fast convergence for a wide selection of expectation values 

TABLE III: The energy contributions to the ground state of 
4 He. These data use values of the physical constants 1/a = 
137.035999679 and m e /m a = 0.000137093355571, where a is 
the fine-structure constant, m e is the mass of the electron, 
and m a is the mass of an alpha particle [4j], y2] . The errors 
do not include the uncertainties in these values. 



Energy 


This Work" 


Drake [35] 


(H ) 


-2.9037243764(8) 


-2.9037243770341195 


\-fJmass/ 




-7.2006570459(3) x 10" 4 


{Hoo} 


-7.4069807(1) x 10~ 6 


-7.40698061439(5) x 10~ 6 


(Hu) 


5.879572027(5) x 10" 4 


5.8795720265(4) x 10" 4 


(Hssc) 


3.55818982(1) x 10~ 5 


3.558189840(7) x 10~ 5 


{Hmp} 


2.18103579(2) x 10~ 5 


2.1810357753732 x 10~ 5 



"Values come from the n = 23 calculation. The errors are calcu- 
lated by assuming an uncertainty five times greater than the fits 
given in Tab. [IT] to account for the spread about these fits. 

''Direct evaluation of the operators pf (i = 1, 2) on the ket yields 
delta function contributions which are unsuitable for direct numer- 
ical evaluation on the grid. So Eq. IC8I cannot be used to produce 
an exponentially accurate expectation value. As is well known, in- 
stead applying p 2 to both the bra and ket produces well-behaved 
functions, but we do not carry out this calculation in this article. 



and matrix elements like the oscillator strength. Varia- 
tional approaches minimize energy-weighted errors but 
generally do not yield comparable results for other op- 
erators. In contrast, we found that the pseudospectral 
method produced errors and convergence rates that were 
very similar for all the quantities studied including en- 
ergy. 

The approach should be widely applicable. No fine 
tuning was done to improve convergence other than en- 
suring non-analytic behavior was treated properly. The 
numerical method we developed was capable of solving 
the large matrix problems with modest computational 
resources. The calculations were pushed to the limits 
of double precision arithmetic. Higher precision floating 
point arithmetic will be necessary to go further. 

This work generalized our previous treatment from S 
to P states and demonstrated the calculation of a variety 
of matrix elements. It can be further extended to higher 
angular momenta in a straightforward manner, albeit at 
larger computational cost. 

The oscillator strength of the helium 1 4 S — > 2 4 P tran- 
sition was calculated to about the same accuracy as the 
most accurate value in the literature [35( and was found 
to agree to the expected precision. 



Appendix A: Bhatia and Temkin Hamiltonian 

Bhatia and Temkin J27[ derived and we checked the 
following explicit expressions that make up the Hamilto- 
nian in their three-three splitting: 
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He 



V 



H 1 



H 1 



K 



lA 1 / d o d 1 d . d 

Z^T ' ~ r > — + ■ /■ */■ sin fli 2 -^- ) 
i=i r * 

Z Z 1 



2 ^ r 2 \dri ,l dn ' sin 9 12 d6 12 



'dOu 



ri r 2 r 12 

cot 0i2 if v = 7 
, (-1)" if ^7 

^stX"' + k2 sin 0i2 - 7 cot 0i2 '( z + ^ if v = 7 

" 1 i/K(2cos0i2+4sin0i2g|-)-?(Z + l)(5i K if !/ 7^ 7 



(1 - s 0K - S lK + (-iy5 2K )hZB lK ,_ 1 
' 2 l - 



H2 Kl = (l-v6 OK )h2B lM 



cot 0i2 if v = 7 
(-l)-r if ^7 



1 / 1 1/7 



>smfi2 \r| r 2 



(1 + <5 2k (V2 - 1))VC - « + 1)(« - « + 2)(Z + k)(Z + k - 1). 

I 



(Al) 
(A2) 
(A3) 

(A4) 

(A5) 

(A6) 
(A7) 



Appendix B: Matrix methods 
1. Formalism 

To solve for the wave function with given k, I and s 
and any m, one must calculate the values of g^ ls for each 
k and v that enters the summation in Eq. 1171 In this 
section we suppress writing k, I, s and m indices; only v 
and k will appear explicitly. There are two types of con- 
ditions which must be satisfied: the Schrodinger equation 
and the boundary conditions. 

The k values of interest are K m , the minimum value, 
K m + 2, ... up to km, the maximum value. The minimum 
and maximum values depend upon parity, I and v (for 
notational clarity omitted). The minimum k is 



v and k. Assemble these in a column vector form that 
enumerates the full set of k for a fixed v 



(B3) 



The length of this column vector is I = 1 + (km — K m )/2, 
which takes on the values \l/2\ or [7/2] . The size of the 
matrix problem increases linearly with I. 

The Schrodinger equation can be represented in matrix 
form: 



/ *L 


\ 


9n m +2 




\ 9km 


) 



y + ^{i-{-i) v k) 



and the maximum is 
Km = 2 



-(1-fc). 



(Bl) 



(B2) 



Let g v K stand for all the grid point values for a given 



J 



H°o 



(H s 
Hi 




= 0, 



(B4) 



where E is the energy, Hs is the S-wave part and iJJ the 
non-S-wave part of the Hamiltonian, and 1 is the identity 
matrix. iJJ and 1 are square matrices with dimensions 
I x I. Explicitly, if J is the tridiagonal matrix 



IP 



( H 1 


H 1 
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H 1 


H 1 


H 1 











H 1 


H 1 









V 







H 1 


H 1 

-1 Hu,K M fl / 



(B5) 



15 



The third subscript on the fl"J K n labels the coupling of 
the individual g functions in k. For the S and P states 
calculated in this article, H'J is only a one by one matrix. 

The pseudospectral matrices H s and HJ Kn (for spe- 
cific v, «, 7 and n) are constructed from Eq. [6] with H 
replaced by Hs or -ffJ Kra , respectively (see appendix [Alfor 
explicit forms of these operators). These single elements 
are large matrices having dimensions set by the number 
of grid points. For multiple subdomains, they are block 
diagonal. The pseudospectral matrix is constructed for 
the subdomain's grid points. The number of columns 
and rows of an element equals the total number of grid 
points in all the subdomains. 

The boundary conditions can be written as 




where 



B„ 



So 
Si 



/ B j» 

Bi m+2 

V 



= 0, 



(B6) 



\ 



B'i" J 



(B7) 



is a diagonal matrix of the same size as i?J, and j m = 
v + K m + 1 + 8 and jm = v + km + 1 + s. Each Bl is a rect- 
angular matrix of the same width as H^ Kn , but a smaller 
height corresponding to the number of grid points near 
internal boundaries or where a symmetry condition holds. 
If j is even (odd) Bl enforces zero derivative (value) along 
the symmetryplane. 

As in Ref. [l| each of the Bl matrices can be split into 
two sub-matrices 



(B8) 



Bl = [Bl x Bl v , 

and similarly splitting the vector g v K 



9 K 




yields the equation 



Hi9* + B J v2 g» 2 = 0, 



(B9) 



(BIO) 



where the vector and matrix have been ordered so that 
the index 1 refers to the rib boundary points and the 
index 2 refers to the rii interior points. The grid point 
nearest to the boundary, at which an explicit boundary 
condition is given is considered a boundary point. B 3 yl is 
an rib by rib matrix and B 3 v2 is an rib by n, matrix. The 
total number of grid points is n t = rib + rii. 

Each n t by rit block of the Hamiltonian matrix iJJ 
(Eq. IB5[) can be split in a similar way, 



H 1 

VKTl 




H 1 H 1 

H 1 H 1 



(Bll) 



TABLE IV: The matrix sizes n m and number of non-zero 
elements unz for each resolution n. 



Resolution 


i 


S States 


i 


P States 


n 


rim 


riNZ 


rim. 


riNZ 


7 


1 512 


182 952 


3 024 


573 720 


8 


2 352 


381 024 


4 704 


1 204 448 


9 


3 456 


722 304 


6 912 


2 297 276 


10 


4 860 


1 273 320 


9 720 


4 069 800 


11 


6 600 


2 118 600 


13 200 


6 798 440 


12 


8 712 


3 362 832 


17 424 


10 826 640 


13 


11 232 


5 133 024 


22 464 


16 571 568 


14 


14 196 


7 580 664 


28 392 


24 531 416 


15 


17 640 


10 883 880 


35 280 


35 292 600 


16 


21 600 


15 249 600 


43 200 


49 536 960 


17 


26 112 


20 915 712 


52 224 


68 048 960 


18 


31 212 


28 153 224 


62 424 


91 722 888 


19 


36 936 


37 268 424 


73 872 


121 570 056 


20 


43 320 


48 605 040 


86 640 


158 726 000 


21 


50 400 


62 546 400 






22 


58 212 


79 517 592 






23 


66 792 


99 987 624 







There are n t + rib equations and n t unknowns (g\ and 
g 2 ) as well as the eigenvalue. One could approximately 
solve these equations with singular value decomposition 
[24| . but it is much faster to simply discard the first rib 
rows of each i?J K „ (one should still check after finding a 
solution that it approximately satisfies those rows of the 
matrix equation) and incorporate the boundary condi- 
tions into the remaining eigenvalue problem by replacing 
each H2 Kn with 



fry _> h"t _ u 1 (B j "r 1 /^ 

n vK,n ^ T1 UK,n22 rL VKn2\\ D v\) D u2' 



(B12) 



where B J vl has an inverse because all of its rows are lin- 
early independent (otherwise more than one boundary 
condition would have been specified for a given bound- 
ary point). Calculating the inverse is computationally 
inexpensive since rib <ti n t . The eigenvector gives g v K2 
and one solves for g v Kl with 

<& = -(Bii)- 1 Bi 2 9:2- (B13) 



2. Matrix Eigenvalue Solution 

The number of grid points in each sub-domain, 
{x, 0,6*12} or {a;,C,/3i2}, was n t = 2n x n x n; greater 
resolution is needed along the semi-infinite coordinate. 
This leads to a Hamiltonian matrix size of rit x nt for S 
states and 2n t x 2rit for odd parity P states. After solv- 
ing for boundary conditions with the above procedure, 



n b +rii 



these are reduced to n r , 
tivcly, where n m — rii = 



x n,, 
3 



(in 



^ and 2n m 
12n 2 + 6n. 



x 2n m , respec- 
Thc number of 
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FIG. 11: (Color online). A log-log plot of the spectral condi- 
tion number c of the pseudospectral matrices as a function of 
resolution. 



non-zero elements n^z scales as n 4 . For n — 20, this cor- 
responds to 560 MB and 1.8 GB, respectively, of memory 
required to store the matrix. 5 The sizes of the matrices 
and the number of non-zero elements is given in Tab. IIVI 

The method of inverse iteration [24[ was used to find 
eigenvalues with a shift equal to the known eigenvalues 
plus 10~ 4 so that the matrix is not too singular. Each 
iteration requires a matrix solve. For the smaller matri- 
ces (up to 17, 000 x 17, 000), these solves were performed 
using Mathematica's [91( multifrontal matrix solve rou- 
tine. This method is fast (eigenvalues can be calculated 
in about 10 minutes for that size) but 8 GB of RAM 
was insufficient to do larger sizes. For larger matrices, 
the generalized minimal residual (GMRES) method of 
PETSc J92|-i3 was used. The GMRES method produces 
a solution with the Krylov space of the matrix and is 
more memory efficient. 

Preconditioning is essential for solving large matrix 
problems. A measure of how hard a matrix problem is 
to solve (how fast a method converges) is the spectral 
condition number, defined as 



I ^max | 
Amin 



(B14) 



where A max and A m j n are the eigenvalues with the largest 
and smallest magnitudes, respectively. The spectral con- 
dition numbers of pseudospectral matrices grow rather 
fast with increasing resolution [95|, [96| . For the problem 
at hand, it is plotted versus resolution in Fig. QTJ It 
starts out large and grows asymptotically as n 12 . An 
ill-posed problem has a condition number which grows 
exponentially [97] . This problem is well-posed but in or- 
der to solve this system of equations preconditioning is 



necessary. A reasonable preconditioner is a matrix pro- 
duced by a second order finite differencing scheme on the 
same set of grid points [23, [95|, [96| . The precondition- 
ing matrix solves are further preconditioned with a block 
Jacobi preconditioner. 

The modified Gramm-Schmidt procedure was used to 
orthogonalize the Krylov subspacc. Furthermore, the 
GMRES restart parameter, m, needs to be very large for 
convergence, empirically, m = 1.3n„i , where n m x n m 
is the matrix size. The computation time scales as n^, 
which for the largest matrix size was about a day run- 
ning on six 2 GHz processors. The eigenvalue solver is 
the slowest part of the entire computation. 

All calculations were done with double precision arith- 
metic. This gives some minimum error in the calculated 
cigenstatc. The effect is relatively big for the small ex- 
ponential tail. The key observation is that the wave 
function no longer decreases at the theoretically expected 
asymptotic rate when it drops to about 10 -8 ' 7 of its max- 
imum value, after which it takes on a seemingly random 
value less than this magnitude. This value is independent 
of resolution because of the limits of machine precision 
arithmetic. It is possible that the asymptotic tail could 
be better calculated with a better preconditioner. 

The issue of the asymptotic behavior is important. 
Since a constant value for the wave function on a semi- 
infinite domain leads to divergent matrix elements, 6 we 
set any value of the eigenvector below this threshold to 
zero. 



3. Quadrature 

In this article, it is necessary to calculate matrix el- 
ements of the form (i|0|j), where \i) and \j) are two 
quantum states and O is some operator. This calculation 
requires numerical integration. Pseudospectral methods, 
by design, use quadrature points as the grid points. A 
one dimensional function /[X] can be numerically inte- 
grated from X=— ltoX = l with weight function g[X] 
by 

J /[X]. 9 [X]dX«^>JpC], (B15) 

where Wi is the quadrature weight specific to the weight- 
ing function g at grid point X 1 . This quadrature formula 
is exponentially accurate with increasing resolution if / 
is smooth over the domain — 1 < X < 1. The problems 
solved in this article are three-dimensional with three 
overlapping subdomains. A separate quadrature can be 
done in each sub-domain. This is illustrated for domain 



5 Note: some eigenvalue solvers do not require one to store this 
matrix and simply require a function which can calculate the 
matrix times a given vector. 



6 For a finite resolution, the quadrature still leads to a finite result 
with an error enhanced by at most 10 4 for the cases calculated 
in this article. 
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0.20 
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(p/n 



FIG. 12: (Color online). This is the arrangement of grid 
points of the three domains at a constant value of p in cj> an d 
#12 coordinates for n = 10. As in Fig. [TJ the blue circles, 
red crosses, and green pluses belong to domains D\, Di, and 
D-i, respectively. Also shown are the overlap grid points in 
D\ n Da (purple stars) and Di n Dz (brown squares). The 
electron-electron singularity is visible at the lower right hand 
corner (solid disk at <f> — 7r/4, #12 — 0) as well as the line of 
symmetry on the right side (dashed line at (f> = 7r/4 where 
r\ — r 2 ). 



so that —1 < Xi,X2,X3 < 1. Integrals over D\ use three- 
dimensional sums analogous to Eq. IB15I Since the ranges 
are fixed, the order of nesting is immaterial. To satisfy 
the requirement that / is smooth (up to the logarithmic 
singularity at p = 0), choose g = l, 7 which corresponds 
to Legendre quadrature points, which are used for all 
calculations in this article instead of Chebyshev which 
were used in Ref. pj. 

If all the subdomains are non-overlapping, then the 
above scheme is sufficient for all integrals. However, no 
set of non-overlapping subdomains for which / is smooth 
could be found. 8 A method is needed for handling over- 
lapping regions, which the above scheme double counts 
if a quadrature is performed in each sub-domain. For 
these regions, an interpolation was performed to two new 
In x n x n grids spanning the overlap regions, shown in 
Fig[T2l For the pseudospectral method, interpolation is 
done to the same order as the grid size. A quadrature 
can then be done over the overlap regions, which are used 
to correct the overall integration. 

The overlap region is divided into two subdomains 



D\ with coordinates {x, <j), 612} and ranges —1 < x < 1, 
< <t> < 1/2, and < 6 12 < tt. Define 



Xi = x 

X 2 = 40-1 

2^12 



(B16) 
(B17) 



X, = 



-1, 



D 13 = DxnDa 
D 23 = D 2 HD 3 . 



(B19) 
(B20) 



(B18) These subdomains satisfy 



D 13 
D 23 



-1 <x< 1, 
-l<x< 1, 



/>min < 4> < 



><i. 



i, < 012 < 012,max[0] 

arccosf < 9 12 < 6»i 2 , max [</>], 



(B21) 



r 



where </> m i n = tt/4 — 1/2 is determined by £ = 1/2 and 
012 = and 9i2, m ax[4>] = arccos[cos 1 esc 2<j>] is deter- 
mined by C = 1/2. One defines appropriate {Xi, X2, X3}. 
For example, in D\ 3 



Xi 



X 2 = 2 



4>- <l>v 



- 1 



X 3 — 20i2 ima x[0]012 — 1- 



(B22) 
(B23) 
(B24) 



Now one calculates the nested sum with X3 innermost 
since the range of #12 depends upon <p. 

The function values at the points necessary for the 
quadrature {x^ 1 , ft 2 , B 3 ^ 3 } ar e calculated with interpo- 



lation 



J 

(B25) 
where Cj refers to the effective basis defined in Eq. [3] and 
J= {ji, J2,J3}- 

Sometimes / involves a Dirac delta function. In such 
a case, one integrates out the delta function analytically. 
One is left with a two dimensional integral on the surface 
where the argument of the delta function is zero. This 
entails first interpolating to that surface using Eq. IB25I 
One can then proceed normally with a two-dimensional 
quadrature. 



18 



Appendix C: Calculating matrix elements with 
Bhatia and Temkin's radial functions 



1. Oscillator Strength 



In the Bhatia and Temkin three-three splitting [271 ] . 
the matrix elements for an *S — >• *P oscillator strength 
transition are written: 



E 



iCSIDrPm)! 2 



/ dT 9000 ( d D9llQ + d D9iw) 



(CI) 
r\r\ sm6i2dridr2d8i2, D is one of the opera- 
tors found inside the matrix elements of Eqs. [22] and the 
operators d l D are given by 



rfi 



4 



4 



4 



(n +r 2 )cos — 

i \ ■ ® 12 

(n -r 2 )sin — 

(ri + r 2 )(3 + cos 6<i2 ) 
4rir2 cos -^ 

2 \ 9ri 9?'2 
(r 1 +r 2 )sin^f 

T\T2 3d\2 

(?'i - r 2 )(-3 + cos6>i 2 ) 
4rir- 2 sin %^ 

+ sin ^i2 (_d d_ 

2 \ dri dr 2 



(C2) 
(C3) 



(C4) 



"a 



4 



(ri — r 2 ) cos -j 2 - 


9 


rir 2 


9^12 


Z(r 2 + r|) cos %^ 




r^2 




Z(rl-rl) sin ^f 





2 2 



(C5) 

(C6) 
(C7) 



2. Expectation Values 

Similarly, an expectation value for an S state is calcu- 
lated by 



CSIDfS) = J dTg 000 d D g° 000 . 



(C8) 



Most of the operators d° D used for expectation values in 
this article have trivial forms. We write here only the 



two most complicated ones: 



U P1P2 



1 



nr 2 



smt/i2 ri- \-r 2 -^ , „„ 



-T\T2 COS6/i2 



if' 



dr\dr2 



COS 6/i2 



^ 12 2 



1 







d°H 00 



sin 6*i2 96»i2 

„2 



(C9) 



a 



S111W12 X12 



<) 



' dr\ 



■x 2 \- 







<9r 2 / 96*12 



+T\T2Z + 



O 2 



dridr 2 



+ Z- 



96-12 2 



12 



i) 



r\T2 smtyi2 OT12 



where 



and 



r? + r\ 2 - ri ■ r 2 



(CIO) 



(Cll) 



z± = (I±3)cos6»i 2 (cos6'i2- / o 2 / 2r -i 7, 2) + sin 2 6li2. (C12) 

All of these forms must be converted to the appropriate 
coordinates in each subdomain. 



Appendix D: History of Oscillator Calculations 

Tabic [V] summarizes the last half century's theoreti- 
cal studies of the nonrelativistic, electric dipole oscillator 
strength. The prime criterion for inclusion in the Table is 
that a numerical value for the oscillator strength for the 
specific transition 1 X S — > 2 : P be calculated and quoted. 
We do not indicate in this Table other transitions cal- 
culated even though these often constitute the bulk of a 
paper's research results. In broadest terms, the entries 
illustrate progress in achieving higher accuracy for the 
specific transition and/or testing new methods designed 
to yield more extensive sets of bound-bound oscillator 
strengths. 

Many methods appearing in Table [V] arc variational 
and utilize the exact interaction pote ntial of the nonrel- 
ativistic Hamiltonian [33, [35|, l98Ml06j . Variational meth- 
ods are especially useful when electron correlation is im- 
portant and ground state properties are sought. There 
are many strategies for selecting bases and suitable vari- 
ational parameters. This flexibility may become cumber- 
some for the study of highly excited states if lower level 
states must be projected out as a preliminary step (e.g. 
if the trial wave function is not linear in the unknown 
parameters and one seeks to enforce orthogonality of the 
excited state with respect to lower states). Errors in the 
eigenproblem accumulate and higher levels are harder to 
find accurately, even when the wave function is linear in 
the variational parameters. 
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A general conclusion is that some basis choices do a 
better job representing the parts of the wave function 
critical to oscillator strength calculat ions . Configur ation 
interaction (CI) calculations [H HI, EM fl07l . Il08t con- 
verge but suffer from the absence of odd powers of the 
inter-electronic distance [44( . Perimetric 1511 coordinates 
3l |35|j4iM IMITol and Hylleraas [If coordinates 
loOMflOa Ill0| include terms of this sort. Sys- 



tematic variational studies using bases incorporating the 
inter-electronic distance have yielded some of the more 
accurate calculations to date. A Hylleraas expansion is 
used by Drake who determined the oscillator strengths 
to seven decimal digits |35( , the most precise calculations 
thus far, as well as some finite-nuclear-mass and rela- 
tivistic corrections. At this stage further nonrelativistic 
calculations of the oscillator strength are probably less 
important than the inclusion of spin-orbit, mass polar- 
ization and low-order relativistic effects. 

Expansions in terms of orthogonal functions often pro- 
duce basis elements of increasing complexity. Alterna- 
tively, one can use larger numbers of simpler functions. 
One important example is the exponential basis [75l |76| | 
(exponential functions of n, r 2 and ry^), which has the 
great advantage of having an easy to calculate Hamilto- 
nian matrix at the expense of violating cusp con ditions. 
This basis was used by Cann and Thakkar |103| to get 
many different oscillator strengths for S — > P and P — > 
D transitions of helium- like atoms. They got the 1 : S 
— > 2 1 P oscillator strength correct to five decimal places. 

The central field approximation |31| is suitable when 
electrons are nearly uncorrclatcd and exchange effects are 
negligible. The essence of this approximation is twofold: 
(1) the multi-electron wave function is written in terms 
of products of one-electron functions and (2) each elec- 
tron experiences a potential which is a function only of 
its distance to the nucleus. The omission of explicit inter- 
electronic coordinates hinders convergence but greatly 
simplifies the variational problem. Green et al. [99( pro- 
duced tables of S — > P and P — >• S transitions using the 
configuration interaction form for the wave functions. 

There exist many different approximations to repre- 
senting the fully correlated wavef unctions. Multiconfig- 
uration Hartree-Fock recovers some but not all of the 
electron correlation energy and yields improved oscilla- 
tor strengths compared to Hartree-Fock treatments |lll| . 
The coupled cluster expansion (roughly analogous to a 
truncated for m of configuration interaction) also yields 
better results J112J . 

Simplifications are frequently made to generate 
comprehensive but approximate oscillator strength 
databases. With this approach the physical as opposed 
to numerical errors may be difficult to gauge. For a two- 
electron atom the Hamiltonian may be written 



H ss H\+ H 2 



Hi = ^ + Ui[n], 



(Dl) 

(D2) 



electron cloud. If the matrix element is dominated by 
the wave function at large distances one may adopt the 
asymptotic form of the p otent ial in that limit to give the 
Coulomb approximation J113l |. 



U, 



Z - 1 



(D3) 



In this approximation the regularity condition at r = 
no longer applies; one needs an alternate method of de- 
termining the discrete energy eigenvalues. These may 
be borrowed from experimental measurements or other 
theoretical calculations and are referred t o as "hybrid" 
results in Table El Wiese et al. |ll4 Ill5l | used this ap- 
proximation (with exchange effects) to calculate oscilla- 
tor strengths for the elements from hydrogen to calcium, 



Cameron et al. 
and Thcodosiou 



11611 tabulated 95 different transitions, 
117| produced extensive tables with er- 
rors better than 10% based on a more sophisticated form 
jl!8l | of Ui. He calculated the oscillator strength of the 
l'S-> 2 ^P tr ansition to four decimal places. Runge and 
Valance [119J ] developed a similar approach based on the 
atomic Fucs potential for the valence electron, 



U[r) 



Z 

r 



E 

i=o 



BiPi 



(D4) 



where Bi is an adjustable parameter and Pi is the pro- 
jection operator onto a subspace of given angular mo- 
mentum I. Currently, the most compl ete t abulation of 
transitions is given by Wiese and Fuhr [l20j . 

The Tabic includes calculations ba sed on pe rturbation 
theory. Sanders, Scherr, and Knight J12l[|l22| developed 
a \jZ expansion, in which the electron-electron interac- 
tion is the perturbation. Even for Z = 2, calculations 
could be carried out to high enough order that the oscil- 
lator strengths converged to three decimal places for the 
helium 1 1 S — > 2 1 P transition. One merit of this approach 
is that it yields oscillator strength as a function of Z and 
with an improving accuracy as Z and/or excitation levels 
increases. 

Devine and Stewart |l23l . ll24| divided the Hamiltonian 
into two parts 



Ho — Hhf + Hi , 



(D5) 



where Ui accounts for the screening of the nucleus by the 



where Hhf is the Hamiltonian projected into the sub- 
space spanned by solutions of the Hartree-Fock type and 
H\ is the difference between this operator and the full 
nonrelativistic Hamiltonian Hq. The operator Hi was 
treated as a perturbation parameter using wave functions 
derived from the frozen Hartree-Fock core. They derived 
oscillator strengths correct to three decimal places using 
second-order perturbation theory. 

Finally, some results do not attempt to calculate os- 
cillator strengths with greater precision or for larger sets 
of transitions. Anderson and Weinhold J110| calculated 
oscillator strengths and rigorous bounds on those values. 
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The 1 1 S —> 2 X P oscillator strength results derived in 
this paper by the pseudospectral method are not listed 
in Table [V] but match the accuracy of the most accu- 



rate included. The method has not yet been tested on 
transitions involving other states. 



J 



TABLE V: A brief history of theoretical calculations of the non- 
relativistic, electric dipole contribution to the oscillator strength for I^S 
— > 2 X P transition. 



Authors 


Method 


Value 


Notes 


Trefftz et al. 


HF, explicit corr 


0.3113L 0.2719V 


Table 4, wf: 2 orbitals with ri2, v 


EU 






preferred 


Dalgarno and 


Sum rules 


0.239 


Table 1, /'s from earlier calcula- 


Lynn 1957 [H 






tions modified for conformity 


Dalgarno and 


var, Hyll 


0.275 


quoted Low and Stewart in Table 2, 


Stewart 1960 






wfs: 6 parameter S, Z* hydrogenic 


na 






P 


Schiff et al. [9g] 


var, peri coord 


0.276159L 


Table I, extrap 56, 120, 220 term 






0.276164V 


wfs, method D 






0.276149A 








0.276154L 


Table VII, extrap 56, 120, 220 term 






0.276150V 


wfs, method C 






0.27616 


Table IX, ±0.00001, summary. 


Green et al. [99l| 


var, CF with exch and CI 


0.27537L 


Table 1, Slater orbitals, wf: 50 






0.27586V 


terms IS, 42 terms 2P, Z* , hybrid 






0.26908A 




Weiss [lOOl] 


var, Hyll coords 


0.2759L 0.2761V 


Table 2, wf: 53 terms IS, 52 terms 
2P; EFs; hybrid 


Cohen and Kelly 


HF, FC; one valence electron; some 


0.112L 


Table V 


EU 


exch 






Dalgarno and 


HF, ~ Z- 1 


0.373 


Table 3, first order in Z _1 


Parkinson [i"2SJ] 








Chong and Ben- 
ston flOlTl 


var, constrained by off-diagonal hy- 


0.26385 


f from M (0.41620, Table II) and 


pervirial theorem 




calculated energy (0.77459), wf: 7 








terms for IS, 2 terms for 2P, Z* 


Sanders and 


var, Hyll coord, Z~ x 


0.276113L 


Table XVIII, wfs: 100 terms, 9-th 


Scherr [ili^ 




0.276182V 
0.276012A 


order in Z 


Cameron et al. 


HF FC, one valence electron 


0.281L 0.255V 


Table I 


Ea 








Schiff et al. [33| 


var, peri coord 


0.276165V 


Table XIV, wfs: up to 1078 terms 
for S state, 364 for P states; con- 
verged to within number of digits 
quoted 


Devine and Stew- 


HF, FC, Pert 


0.2760L 0.2749V 


Table 2, iterated result, wf: 77 


art EH 




0.2771A 


terms for S state, 65 for P state 


Laughlin EU 


Z^ 1 , mod screening 


0.29834 


Table 4, f from expansion 
coefficients 


Anderson and 


rigorous limits 


0.2747 - 0.2775 


Table IV 


Weinhold EH 








Froese Fischer 


MCHF 


0.2753L 0.2744V 


Table 2 


m 








Leopold and Co- 


upper bounds 


< 0.29678 


bound from a 2 (Table 1) and best 


hen Em 






NR energy; hybrid 


Davis and Chung 


CI, no rl2 corr, AMPW 


0.2721L 0.2758V 


Table V, 110 terms S and P states 











Roginsky and 


Modified wf 


0.256L+V 


Table 1, product wfs with Z* 


Klapisch Ell 








Kono and Hattori 


var, double Hyll, ECFs 


0.27616 


Table III, 138 terms S and 140 


Eol 






terms P, 3 nonlinear parameters (2 
set, 1 optimized), ±0.00001 


Continued on next page 
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TABLE V — continued from previous 


page 


Authors 


Method 


Value 


Notes 


Theodosiou [133J 


Valence electron in potential 


0.2761L 


Table I, HF Slater potential, hybrid 
(experimental) 


Park et al |134| 


HSA 


0.291L 0.342A 


Table I, initial (final) wf 4 (6) an- 
gular momentum pairs 


Theodosiou |117| 


as above 


0.27643L 


Table I 


Fernley et al. 


CC expansion 


0.2811 


Table 3, Is, 2s, 2p, Id, 3p one- 


EH 






electron states and product states; 
R-matrix inner region, numerical 
integration outer region 


Sanders and 


var, Hyll, ~ Z _1 , pert 


0.27774 


Table V, wfs and energies from |l29|] 


Knight [1221] 








Abrashkevich et 


HSAnacc 


0.2763L 0.2844A 


Table 2, initial (final) 6 (4) radial 


al. [13511 






equations, 100 finite elements 


Cann and 


var, exp, ECFs 


0.27617 


Table V, 100 terms, 6 nonlinear pa- 


Thakkar [l03| 






rameters (error of 0.7 — 2.99 units 
in last digit) 


Tang e£ a/. |l36l] 


HSCC CC 


0.2762L 0.2763A 


Table I 


Chen [lOj| 


CI with B-splines 


0.27611 


Table 12, 150 9-th and 10-th order 
splines for S, 137 for P, uncertainty 
< 0.01% 


Chen [1051] 


CI with B-splines 


0.276163L 


Table 13, 150 9-th and 10-th order 






0.276076V 


splines for S, 147 for P, hybrid (best 
NR energies) 


Yang [lOefl 


MELL, peri 


0.276165L 


Table 3.7, 680 terms, 2 nonlinear 






0.276165V 


parameters 


Drake ^ 


var, double Hyll 


0.2761647 


Table 11.11, nonlinear scale 
parameters 


Masili e£ al [l37t| 


HSAnacc 


0.2761957 


Table 4, initial (final) 13 (15) radial 
equations 


Alexander and 


varMC 


0.2761L 0.2706V 


Table V, largest wfs, rotated 


Coldwell [1381] 




0.2758A 


method 



TABLE VI: Abbreviations in above table. 



Symbol 


Meaning 


A 


acceleration form for oscillator strength 




AMPW scale params 


angular momentum partial wave scale parameters 




CC expansion 


close coupling expansion 




CF 


central field (no separation coordinate) 




CI 


configuration interaction 




corr 


correlation factors 




double Hyll 


double Hylleraas coordinate basis 




ECFs 


multiple exponential correlation factors 




EFs 


multiple exponential factors 




exch 


exchange interactions 




extrap 


extrapolation based on 




exp 


exponential basis (exponentials of Hyllerass coordinates) 




f 


oscillator strength 




FC 


frozen core 




HF 


Hartree Fock 




HSA 


adiabatic Hyperspherical coordinate representation 




HSAnacc 


HSA with non-adiabatic channel coupling 




HSCC 


Hyperspherical coordinate representation; CC expansion 




hybrid 


energy not taken from parmeterized wave function; input 
ment or other calculations 


from experi- 


Hyll coord 


Hylleraas coordinates 




L 


length form for oscillator strength 




mod screening 


modified screening approximation 




M 


dipole moment 




MCHF 


multiconfiguration Hartree Fock 






Continued 


on next page 
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TABLE VI — continued from previous page 




Symbol 


Meaning 


MELL 




matrix expansion in exponentials, Laguerre polynomials and 
tions of total orbital angular momentum 


3igcnfunc- 


NR 




non-relativistic 




peri coord 




perimetric coordinates 




Pert 




perturbation theory corrections 




var 




variational 




varMC 




variational Monte Carlo 




V 




velocity form for oscillator strength 




wf, wfs 




wave function, wave functions 




z- 1 




expansion in inverse powers of Z 




-z- 1 




expansion in inverse powers of Z with additional corrections 




z* 




nonlinear, effective nuclear charge parameter 





r 
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